Controlled-source electromagnetic (“CSEM”) surveys are a useful geophysical tool for evaluating the presence of hydrocarbon-bearing strata within the earth. CSEM surveys typically record the electromagnetic (“EM”) signal induced in the earth by a source or transmitter, and measured at one or more receivers. The behavior of this signal as a function of transmitter location, frequency, and separation (or offset) between transmitter and receiver can be used to determine the properties of the rock, specifically whether hydrocarbons are present. One of the properties is electrical resistivity. Thus, CSEM measurements are typically used to determine the spatially varying resistivity of the subsurface region of interest.
CSEM surveys may be performed over land or water. In a water environment, CSEM data are typically acquired by towing an electric dipole transmitter antenna behind a vessel. The seafloor above the region of interest is populated with a plurality of receivers. The receivers typically have multiple sensors designed to record different vector components of the electric and/or magnetic fields. The receivers are placed on the seafloor at distances of two to three kilometers apart by the same vessel that tows the transmitter or by another vessel. The receivers are weighted, and thus sink to the seafloor. After placement, the vessel broadcasts a low frequency dipole signal through the transmitter. Typically, the signal ranges from ⅛ to 2 Hertz, but may be as low as 0.001 Hertz and as high as 10 Hertz. The low frequency allows signal penetration into subsurface regions below the seafloor. After the survey is completed, a command, e.g. an acoustic command, from the tow vessel or another vessel causes the receivers to release their weights, and thus to rise to the surface. The receivers, with their data, are then retrieved. To form the data, the transmitter antenna is typically towed a few tens of meters above the seafloor by the vessel on the sea surface. The antenna sends out a plurality of pulses, and each receiver records the values of the electric and/or magnetic fields in the region around the receiver.
Alternative configurations include stationary transmitters on the seafloor, or in the water column, as well as magnetic transmitter antennae and passive sources, such as magnetotelluric energy. Such energy may be generated, for example, by the interaction of the solar wind with the earth's magnetic field. The transmitting and receiving systems typically operate independently, meaning without any connection between them, thus the receiver data needs to be synchronized with shipboard measurements of transmitter position by comparing clock times on the receivers to time from a shipboard or GPS (Global Positioning System) standard.
The transmitted signal drives electrical current into the subsurface and the water column, resulting in an electromagnetic field throughout the subsurface and the water column, based on the laws of physics encapsulated, for example, in Maxwell's equations. The data collected by each of the receivers is essentially an expression of the amplitude and phase of the resulting electric and/or magnetic fields versus the distance between the receiver and the transmitter. The amplitudes peak when the transmitter is closest to the receiver and are smallest when the transmitter is far away from the receiver. The resulting electric and magnetic fields may vary depending upon whether the subsurface region includes conductive or resistive bodies. Conductive bodies may include, for example, brine-saturated sandstones or shales; whereas resistive bodies may include salt domes, basalt flows, limestones, or hydrocarbon-saturated reservoirs.
CSEM data is typically analyzed in the temporal frequency domain with each signal representing the response of the earth to electromagnetic energy at that temporal frequency. Temporal frequency domain means the data are transformed, typically by Fourier transformation, such that the dependence of the data on time becomes a dependence on frequency. In raw CSEM data, the strength of each frequency component varies depending on how much energy the transmitter broadcasts. That is, the strength depends upon the amplitude and phase of each component in the frequency spectrum of the transmitter, and also on the receiver sensitivity at that frequency. These transmitter and receiver effects are typically removed from the data prior to interpretation.
Note that the resistivity of the subsurface regions, such as the Earth, is generally anisotropic. Thus, the ratio of electric field to applied current depends both on the direction of the current and on the direction of the field. It is further recognized that the subsurface regions may be equivalently characterized by resistivity or its inverse, electrical conductivity. Areas that comprise brine-saturated rock have lower resistivities and higher conductivities than areas that may contain nonconductive materials, such as oil.
Frequency-domain inversion is commonly used to infer resistivity values from CSEM data. See, for example, Carazzone, J. J., et al., “Three Dimensional Imaging of Marine CSEM Data,” Expanded Abstracts, 75th Annual Int'l Meeting, Society of Exploration Geophysicists, pp. 575-578, (2005); and MacGregor, L., et al., “Derisking Exploration Prospects Using Integrated Seismic and Electromagnetic Data—a Falkland Islands Case Study,” The Leading Edge, 26, pp. 356-359, (2007). Inversion is useful because it produces resistivity models of the subsurface consistent with the measured data, specifically, with the amplitude and phase of one or more measured components of the electric and/or magnetic fields at one or more frequencies and receiver locations. Using Maxwell's equations, one-dimensional (1-D), two-dimensional (2-D), and three-dimensional (3-D) inversions of the data generate 1D resistivity profiles, 2D resistivity cross-sections, and 3D resistivity volumes, respectively. In some cases, it may be desirable to obtain and incorporate measurements of the CSEM magnetic fields as well as natural-source electromagnetic data, such as magnetotelluric data, in the inversion together with the CSEM data.
Linearized, iterative inversion of electromagnetic data is discussed, for example, in Commer, M., and Newman, G. A., “New Advances in Three-Dimensional Controlled-Source Electromagnetic Inversion,” Geophys. J. Int., v. 172, pp. 513-535, (2008). This inversion process determines a single resistivity model that best fits the collected data, while simultaneously obeying a smoothness constraint and/or other constraints, such as a priori subsurface resistivity information. Beginning with an initial resistivity model the inversion process forward-simulates electromagnetic data corresponding to the source and receiver configurations in the actual survey. An objective function is then computed, which is a measure of the mismatch between the measured and synthetic data and of the quality of the constraint. If the objective function is sufficiently small, then the inversion process is complete and the resistivity model is accepted as an approximate representation of the actual resistivity of the subsurface region. More commonly, the resistivity model may be adjusted to reduce the objective function. In 2-D and 3-D inversions, the most practical technique is the use of a gradient-based method which determines a model adjustment that is orthogonal (in the space of model parameters) to the gradient of the objective function. (For a discussion of gradient-based methods see Dennis, J. E., Schnabel, R. B. “Numerical Methods for Unconstrained Optimization and Nonlinear Equations,” Printice-Hall, Englewood Cliffs, N.J. (1983) (ISBN-0-13-627216-9). This model adjustment is applied to the resistivity model and the process is repeated, with new synthetic data being generated and a new model adjustment being determined. As the process is repeated (or iterated), the influence of the constraint term is commonly decreased, so that it does not overwhelm the mismatch between measured and synthetic data in the objective function. After a number of iterations, the objective function may either be sufficiently small or no model adjustment can be found that may further reduce the objective function. A resistivity model for which no adjustment can be found that further reduces the objective function is known as a “local minimum” to those practiced in the art of inversion.
An advantage of gradient-based inversion algorithms is that they are relatively computationally inexpensive. However, it is well known that such inversion methods do not typically find the optimal resistivity model. Electromagnetic inversion is a so-called ill-posed problem, which implies that there are in fact a multiplicity of resistivity models that are consistent with the measured electromagnetic data. Therefore, the results of the inversion process are not unique, but depend upon a number of factors, one of which is the initial resistivity model. The inversion process does not in general converge to the global minimum of the objective function but, rather, it converges to a local minimum lying near the location in model space occupied by the initial resistivity model, as discussed with reference to FIG. 1.
FIG. 1 depicts an example where an inversion process settles at a local minimum rather than at the global minimum. The graph 100 depicts a plot of misfit versus model parameter. The graph includes a plurality of local minima, at points 101a-101c, with one of the local minima being the global minimum point 101a. The global minimum point 101a represents the optimal (e.g., best or correct) model with the least data misfit. The vertical axis in the graph 100 represents the data misfit X2, given by equation (1):
                                          X            2                    =                                    1              M                        ⁢                                          ∑                                  i                  =                  1                                N                            ⁢                                                          ⁢                                                                    (                                                                  E                        i                        data                                            -                                              E                        i                        mod                                                              )                                    2                                                  θ                  i                  2                                                                    ,                            (        1        )            where E represents electric and/or magnetic field values and the sum is over frequencies, vector components and source-receiver separations, which are represented by index i. The labels ‘data’ and ‘mod” refer to measured data and modeled data, respectively. The total number of data points is given by M, and the quantity θi is the measurement uncertainty of the ith data point. The equation (1) is a commonly used representation of data misfit. Please note that the measurement of uncertainty is typically represented by the symbol σ; however, we also typically use σ to represent conductivity, so for clarity we use the symbol θ to represent measurement uncertainty.
The inversion process minimizes an objective function that may include the data misfit X2 alone, or the data misfit X2 plus additional regularization term(s) designed to emphasize a certain characteristic of the model, such as its smoothness. In CSEM inversion, the model parameter being solved for is the subsurface conductivity, σ, which is the inverse of the resistivity. For example, in FIG. 1, the horizontal axis represents the value of the conductivity. For each frequency, electromagnetic field component, and source-receiver configuration, Eimod is determined from the model comprised of the totality of the σ values by solving Maxwell's equations.
In the space of model parameters, gradient-based inversion methods determine a path that decreases the data misfit X2. For example, given an initial model estimate 102, gradient inversion moves the model parameter along the path 106 toward the minimum 101c at which point the inversion process halts. As used herein, the “conjugate gradient method” is a method for finding the nearest local minimum (e.g., the smallest value of a set, function, etc., within some local neighborhood) of a function of n variables which presupposes that the gradient of the function can be computed. A gradient in 3-D is, for example, the vector sum of the partial derivatives with respect to the three-coordinate variables x, y, and z of a scalar quantity whose value varies from point to point. In resistivity inversion, the relevant gradient is the n-dimensional vector of derivatives of the objective function with respect to the model conductivities. The conjugate gradient method uses conjugate directions instead of the local gradient for approaching the local minimum. If the vicinity of the minimum has the shape of a long, narrow valley, the minimum is reached in far fewer steps than other methods, such as the method of steepest descent. Because the initial point 102 lies near a local minimum 101c and not the global minimum 101a, the process finds a solution at the local minimum 101c, and not at the global minimum 101a. However, if the initial estimate were chosen to be point 104, then the inversion process finds a solution at the global minimum 101a by traveling along path 107 during the iterations.
Several techniques are available for overcoming the problem of finding local minima instead of the more desirable global minimum. Collectively, these are known as “stochastic” inversions because they attempt to avoid becoming trapped in local minima by employing some randomness in their search patterns. A stochastic process is one whose behavior is non-deterministic in that a subsequent state is determined by the process's predictable actions and by a random element. Thus, stochastic inversion methods may become untrapped at local minimum 101b and may converge to the global minimum at 101a. One stochastic technique is crude Monte-Carlo forward modeling, which involves searching the entire model space, and thus may involve too much computer time for a 3-D problem.
Non-limiting examples of other stochastic techniques include genetic algorithms, simulated annealing, and Markov-Chain Monte-Carlo (MCMC) (which is a method utilizing importance sampling), which may be faster techniques for attempting to find the global minimum. First, the genetic algorithms involve either randomly changing the set of inversion parameters found at each iteration or merging the set with sets of parameters from another iteration to form the parameter set or genetic code for the next iteration. The randomization of the parameters enables the method to avoid local minima Second, the simulated annealing methods mimic the physical process of cooling. They typically start by allowing the parameter set to vary relatively widely in the initial iterations, and then gradually decrease the allowable parameter space as iterations continue. Thus, the parameters are not likely to settle onto a local minimum in the first stages of the iteration, and should eventually converge to the global minimum as the iterations proceed and the allowable parameter range decreases. Finally, the Markov-Chain Monte-Carlo methods combine a random search with a prescription for ensuring that the search explores the correct portion of parameter space. These methods are further described in Chen, J., and Dickens, T., “Effects of Uncertainty in Rock-Physics Models on Reservoir Parameter Estimation Using Marine Seismic AVA and CSEM Data,” Abstracts of the 77th Annual Int'l Meeting, Society of Exploration Geophysicists, pp. 457-461, (2007); Chen, J., and Dickens, T., “Effects of Uncertainty in Rock-Physics Models on Reservoir Parameter Estimation Using Seismic Amplitude Variation with Angle and Controlled-source Electromagnetic Data,” Geophysical Prospecting, 57, pp. 61-74, (2009); and Chen, J., et al., “A Bayesian Model for Gas Saturation Estimation Using Marine Seismic AVA and CSEM Data,” Geophysics, 72, pp. WA85-WA95, (2007). The prototype of these methods is the Metropolis-Hastings algorithm, which makes a random update to some model parameter at each step, and then accepts or rejects the update according to a specified probability distribution. All updates that decrease the misfit are accepted, and a portion of the updates that increase the misfit are allowed, giving the method the ability to jump out of local minima. These are a few examples of stochastic techniques.
Regardless of the method, even with their increased efficiency, these stochastic techniques may be limited because they involve unacceptable computer time to be applicable to most problems in 3-D and many problems in 2-D. That is, a characteristic common to stochastic techniques is that they may require several thousand or even tens of thousands of forward data simulations and evaluations of the objective function for different underlying earth models. A related characteristic of stochastic methods is that they result in a distribution of possible earth models and a measure of the likelihood, or probability distribution, of each. As such, these techniques are computationally expensive.
An additional approach is to utilize a hybrid method. Hybrid methods combine stochastic techniques with deterministic methods. These methods typically interleave iterations of deterministic inversion with stochastic inversion in some fashion designed to restrict the randomness of the stochastic search by exploiting the guaranteed decrease in the objective function from deterministic methods. However, similar to the stochastic techniques, these inversion methods are typically computationally expensive compared to direct inversion algorithms, which attempt to find a minimum in parameter space in as few steps as possible.
Direct, deterministic inversion algorithms include the simplex, steepest descent, conjugate gradient, quasi-Newton, Gauss-Newton, and full Newton techniques. Each of these algorithms performs only a few modeling steps to determine the location of the next parameter set in model space. Typically, these algorithms are designed to find the inversion result using the smallest possible number of forward modeling steps. Each algorithm bases its iterations on local information, so each is guaranteed only to find a local minimum in the solution space. In contrast to stochastic methods, gradient techniques typically only use a few tens or perhaps up to one hundred forward simulations from different underlying earth models. Note that typically the high cost of forward modeling in 2-D and 3-D, plus the increase in the number of parameters required to describe the model, cause stochastic inversion to be infeasible for such problems.
The degree to which 2-D and 3-D problems are underdetermined, and therefore the number of incorrect local minima, is even larger than for the 1-D case. Most areas of interest for hydrocarbon exploration are of such complexity that a 1-D model does not adequately describe them. Accordingly, 2-D and primarily 3-D models are utilized to properly describe the area of interest. Because of the computational requirements mentioned above, inversion of field data is therefore limited to gradient-based methods, so the results depend significantly on which local minimum lies near the starting model.
With presently available algorithms and computers, the application of stochastic inversion techniques to find the optimal global minimum for such problems is not practical. This severely limits the usefulness of inversion for making business decisions. If the inversion produces a result corresponding to a local minimum lying far from the global minimum, the predicted rock properties and/or structural framework may be so inaccurate as to cause the analyst to rely upon or infer inaccurate interpretations from the data. Regardless, even if the global minimum is found, its validity may be called into question, as it does not provably lie at the global minimum.
As stated earlier, one aspect of inversion techniques is that they involve an initial estimate of the earth resistivity model (e.g., a starting model or initial model). The creation of a good starting model is particularly challenging for 3-D models, because large model data sets are manipulated and populated with anisotropic resistivities, while maintaining enough consistency with the large volume of measured data and with other geophysical data to minimize errors and artifacts in the final inverted result. The inversion result is highly dependent upon the starting model used. As such, the 3-D CSEM inversion problem is typically ill-posed because different starting models may lead to different results. Furthermore, the use of a non-optimal starting model can lead to the presence of nonphysical artifacts whose sources may be difficult to identify.
One widely used technique for producing starting models for inversion processes is to make a simplifying assumption that allows the analyst to create a reasonable model. For example, velocity or resistivity models are often initialized by using information from a well, which is then extrapolated through the remainder of the model assuming that velocities or resistivities follow the same compaction trend starting at the water bottom. Similarly, anisotropy parameters may be extrapolated along the structural boundaries of a shale package. Thus, simplified models based on realistic geological assumptions may be useful to cure the deficiencies of nonoptimal starting models.
In other cases, a local layer cake or one and half-dimensional (1.5-D) model may be assumed at each data collection location. The data is analyzed at each location to determine the 1.5-D model that optimally fits the data. These 1.5-D models are then interpolated/extrapolated to produce a 3-D model that incorporates the information available near data-collection locations. For example, see Wahrmund, Leslie A., et al., “Method for Obtaining Resistivity from Controlled Source Electromagnetic Data,” U.S. patent application Ser. No. 12/280,330 filed Aug. 21, 2008 from PCT Application PCT/US2007/004111, Published on Nov. 8, 2007. The terms “1.5-D” and “2.5-D” are commonly used by those familiar with geophysical and processing and inversion methods, and refer to earth models that, while encompassing three spatial dimensions, are described by properties that vary only in depth, or in depth and in one lateral dimension, respectively.
In addition to resistivity information obtained by inverting the electromagnetic data, other types of data may be used to help constrain the resistivity model. Well log data may be used to populate the model near a given well. If seismic reflection data are available for the location, the seismically identified reflectors can be used to correlate the resistivity data between well locations. Information on lithology and velocity gleaned from the seismic data can be used to loosely constrain the initial resistivity values. Again, if there are both well resistivity data and seismic velocities, these can be correlated and used to populate the resistivity model away from the well.
Alternatively, it is possible to use higher-dimensional information, such as a 2-D or 3-D seismic image and a rock-physics model that connects the seismic and electrical properties of the subsurface to initialize an approximate resistivity model in higher dimensions. See Chen, J., and Dickens, T. (2009). The rock-physics model may itself be calibrated by well log measurements or even replaced by a heuristic model based on well-log data.
The foregoing discussion of need in the art is intended to be representative rather than exhaustive. A technology addressing one or more such needs, or some other related shortcoming in the field, may be beneficial in developing geological models representative of subsurface regions and may facilitate prediction of a subsurface location for hydrocarbons based at least in part on such models.
Other related material may be found in the following publications: Eide et al., “Stochastic Reservoir Characterization Conditioned on Seismic Data,” Geophysics, vol. 69, No. 4, pp. 824-836; Luther White et al., “Stochastic Fluid Modulus Inversion,” Geophysics, vol. 67, No. 6, pp. 1835-1843, (November-December 2002); Debeye et al., “Stochastic Inversion”, The Strategic Importance of Oil And Gas Technology Proceedings of the 5th European Union Hydrocarbon Syposium, Edinburgh, U.K., v. 1, pp. 166-175 (1997); U.S. Pat. Nos. 6,067,340, 6,278,948 and 6,970,397; and U.S. Patent App. Pub. No. 2009/0006053.